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The transition interface sampling (TIS) technique allows to overcome large free energy barriers 
within reasonable simulation time, which is impossible for straightforward molecular dynamics. 
Still, the method does not impose an artificial driving force, but it surmounts the timescale problem 
by an importance sampling of true dynamical pathways. Recently, it was shown that the efficiency of 
TIS to calculate reaction rates is less sensitive to the choice of reaction coordinate than those of the 
standard free energy based techniques. This could be an important advantage in complex systems 
for which a good reaction coordinate is usually very difficult to find. We explain the principles of 
this method and discuss some of the promising applications related to zeolite formation. 



I. INTRODUCTION 



Gaining insight in the zeolite formation has not only 
fundamental scientific importance, but could also accel- 
erate momentous technological developments. The appli- 
cations of zeolites are uncountable ranging from cracking 
catalysis of crude oil, gas separation, detergent builders, 
and sensors for pharmaceutical formulations. The spe- 
cific catalytic properties of zeolites lie in their unique 
open crystalline structure that incorporates cages or 
channels with typically nanoscale diameters. The growth 
of the open structure silicon dioxide polymorphs is me- 
diated by so-called template molecules that can be re- 
moved out of the zeolite pores after the crystallization 
process. Besides template molecules, solvent, Si/Al ra- 
tio, temperature, pH, and many other factors play a role 
in determining which zeolite topology is formed. As each 
structure and composition has its unique catalytic prop- 
erties, the synthesis of new zeolite materials has been 
an important branch of chemical research. This devel- 
opment has progressed mainly on the basis of trial-and- 
crror and 'chemical intuition' as a fundamental under- 
standing of zeolite formation is lacking. The clear so- 
lution synthesis studies of silicalite-1 zeolites initiated 
by Schoeman et ai^ were an important step forward for 
the experimental analysis. The use of clear solutions in- 
stead of gels made the analysis of zeolite synthesis much 
more accessible by experimental techniques. Since then, 
this model system has been subject of many studies in- 
cluding x-ray and neutron scattering, infrared (IR) spec- 
troscopy, Nuclear magnetic resonance (NMR), and dy- 
namic light scattering (DLS). These studies revealed that 
upon mixing tetraethylorthosilicate (TEOS), tetrapropy- 
lammonium hydroxide (TPAOH) and water at a cer- 
tain ratio at room temperature sub-colloidal particles are 
formed of several nanometers. Using the freeze drying 
technique^, these particles have been extracted from the 
solution and examined by various techniques such as solid 
state NMR, Fourier Transform IR (FTIR), transmission 
electron microscopy (TEM), and atomic force microscopy 
(AFM). Various models for the structure of these parti- 
cles have been proposed ranging from amorphous bod- 



ies^ to precise framework fragments^. 

The formation of crystalline zeolite particles is initi- 
ated when this suspension is heated upto temperatures 
of 350 K. Light scattering experiments show that the in- 
tensity scattered by the suspension increases only slowly 
in time during the first period of the synthesis. This is 
then followed by a sharp increase, indicating the starting 
point of growth of what will become the final crystals. 
The first period can be associated to a nucleation pro- 
cess, in which a particle has to be formed with a size 
beyond its critical nucleation radius. The formation of 
zeolites consist hence of several stages. First a polymer- 
ization process which eventually leads to the formation 
of sub-colloidal particles, second the nucleation process, 
and finally the crystal growth. 

One of the difficulties in the investigation of the ze- 
olite formation process is that the relevant lengthscales 
of the zeolite formation lie just in between the accessible 
lengthscale of NMR and diffraction techniques''. More- 
over, it is unclear if freezed-dried extractions are identical 
to the silicate particles existing in solution. Since many 
experiments do not allow unequivocal interpretation, it 
is not a surprise that several crystallization mechanisms 
have been proposed. These theories concentrate on the 
structure and shape of the colloidal particles, how these 
particles are formed and how these particles finally con- 
tribute in the formation of the zeolite crystal. 

An example is the nanoslab hypothesis that was pos- 
tulated by some of us. It was inspired by several ex- 
perimental observations^^. This theory assumes that 
at an initial stage precursor particles are formed that 
consist of 30 to 33 Si atoms enclosing a single template 
molecule. These precursors stick together in a block 
shaped particle, the nanoslab, that has already the cor- 
rect crystalline structure. These particles finally form 
the zeolite by a 'clicking-mechanism' when the solution 
is heated up. Others have claimed that the apparent 
evidence of the Si-30/33 precursor particle should be 
attributed to other Si containing species^ or that the 
nanoshaped particle is actually an amorphous identity 
with a layer of template molecules around it4£. Also 
the role of the nanosized particles for the nucleation and 
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crystal growth has been subject of debate. According 
to some groups, the nanoparticles add one by one to 
the growing crystal o 11 ' 12 . Others regard the particles as 
monomer reservoirs: monomer dissolves into the solution 
and attaches to the growing nucleu a 13 ' 14 . Recent publi- 
cation o 14 ' 15 state that an aggregative growth mechanism 
of discrete nanoparticles may dominate the early stage 
of the growth process. However, after a certain size is 
attained, the growth mechanism seems to switch to ad- 
dition of low molecular species, probably monomers. 

In conclusion, despite many years of abundant exper- 
imental research, zeolite synthesis still contains many 
mysteries. Therefore, this field of research is a prototype 
example where computer simulations could give invalu- 
able information. However, before truly realistic simu- 
lations of all stages in the zeolite synthesis can be per- 
formed, a long way has to be gone. Reason for the dif- 
ficulty is that the typical system sizes and timescales at 
which the zeolite formation takes place are generally be- 
yond the capabilities of present computer resources. For 
a correct modeling of the nucleation process, the sim- 
ulation box should at least be larger than the critical 
nucleus. A requirement that is out of limits for quantum 
mechanical calculations and demands the development of 
accurate reactive forcefields. 

So far fully quantum mechanical calculations using 
Density Functional Theory (DFT) have been applied to 
silica polymerization cluster a 16 i 17 . These studies showed 
a stronger stability of silicate 6 rings and linear poly- 
mers compared to smaller rin gs a nd branched polymers. 
Ph effects were considered in jl8l.!9] by analyzing nega- 
tively charged silica clusters that are favorable to neutral 
ones in an alkaline environment. This study revealed 
that internal cyclization if preferred over further linear 
growth^. Barriers for oligomerization were significantly 
reduced for single charged cluster compared to neutral 
one*!*. 

Based on ab initio calculations or experimental 
data, several classical forcefields have been devel- 
ope d 20 ' 21 ' 22 ' 23 ' 24 . These potentials allow the study of 
larger systems including solvent and template molecules. 
However, the existing potentials are not yet very accu- 
rately describing the breaking and making of chemical 
bonds, which presumably requires complex many-body 
terms and polarizable forcefields. Still, studies using 
these approximate potentials can give valuable insights. 
For instance, classical molecular dynamics (MD) simu- 
lations have shed some light to the role of solvent and 
template molecule o 25 ' 26 . These simulations showed that, 
contrary to fully formed cages and rings, open structures 
collapse in the presence of solvent, unless it contained 
strongly bonded template molecules. The early stages 
of silica polymerization dynamics were studied by Rao 
and Gelb^l at high temperatures > 1500 K. These al- 
leviated temperatures were required as upto 600 K, no 
polymerization reaction could be observed within the 
nanoseconds simulation periods. They found that both 
the monomer incorporation and the cluster-cluster aggre- 



gation were important mechanisms for diluted solutions, 
while the first mechanism was dominant in the concen- 
trated systems. Using a implicit solvent model not in- 
cluding template molecules, Wu and Deem analyzed the 
free energy barriers and critical cluster sizes as function 
of pH and Si-monomer concentration at ambient condi- 
tions using a series of advanced Monte Carlo (MC) tech- 
niques 2 ^. They found that the critical clusters for the 
polymerization contained relatively few (« 30 — 40) Si 
atoms. No attempt was made to derive reaction rates 
by calculating transmission coefficients. Even larger sys- 
tems and timescales have been simulated using lattice 
models 2 ^ and kinetic MC (KMC) 3 ^. Relative rates for 
different crystal growth mechanisms via kink and edge 
sites can be derived by mimicking atomic force micro- 
graphs via atomistic simulations 32 ' 33 ' 34 . Still, even KMC 
simulations are usually restricted to growt h 35 ' 36 . The 
time, before the critical nucleus of a zeolite is formed, 
is still too long even for this ultrafast type of dynamical 
simulations. 



It is clear that the simulation methods have made sig- 
nificant progress in recent years. At the early stages of 
Si polymerization, fully quantum mechanical MD stud- 
ies are in our reach using Born-Oppenheimer— or Car- 
Parrinell o 37 ' 38 simulations. Thanks to newly developed 
potentials and coarse grained systems, simulations ap- 
proach the system sizes that are needed to describe the 
template directed zeolite synthesis in solution. Nonethe- 
less, each stage in the zeolite synthesis involves signif- 
icant reaction barriers. This makes the chance to ob- 
serve important reactive events at experimental condi- 
tions within the duration of the simulation period highly 
unlikely The reaction itself is usually very fast and could 
fit perfectly within the window of timescales that are at- 
tainable by the simulation method. However, the system 
will likely spend extensively long periods within the well 
of the reactant state without any reactive event taking 
place. It is, therefore, important to have a method that 
focuses the costly simulation time on the important but 
rare reactive events, while limiting the superfluent explo- 
ration of the reactant well. In this article, we review such 
a method, the transition interface sampling (TIS ) 39 ' 40 ' 41 
method, that allows to concentrate only on those trajec- 
tories that are important for the chemical process. More- 
over, the TIS technique can calculate the frequency for 
occurrence of these successful trajectories within an in- 
finitely long straightforward simulation. Hence, TIS al- 
lows the determination of the rate of the rare event. 



The aim of this article is not to give a fully detailed 
theoretical derivation of the method. This has already 
been published elsewher o 39 ' 40 ' 41 . The goal of this article 
is to give an educative overview of the practical algo- 
rithms and their possible applications related to zeolite 
studies rather than on mathematical aspects. 
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II. TRANSITION INTERFACE SAMPLING 
A. historic perspectives of rare event simulations 

The first theories for treating rare events from a 
microscopic perspective where pioneered by Eyring 4 ^, 
Wigner— , and Horiutu 4 ^ about 20 years before the first 
MD simulation was performed 45 . They introduced the 
concept of Transition State (TS) and the so-called TS 
Theory (TST) approximation. Later Keck demonstrated 
how the TST approximation can be made exact by a dy- 
namical correction, the transmission coefficient 46 . The 
actual application of these theories for molecular simula- 
tion was directed by the works of Bennett^ and Chan- 
dler—, which have made this a standard approach in 
molecular simulation. A crucial point in this reactive 
flux (RF) method is the definition of a suitable reaction 
coordinate (RC). As a first step, the free energy needs 
to be determined along this RC using importance sam- 
pling techniques such as Umbrella Sampling (US) 4 ^ or 
Thermodynamic Integration (TI)^. This result alone is 
sufficient to obtain the TST approximation of the rate, 
which is an upper limit for the actual rate. In the second 
step, the correction to this approximation can be calcu- 
lated by releasing dynamical trajectories from the top of 
the free energy barrier. Only when both steps are com- 
pleted, the exact reaction rate can be calculated. Both 
the free energy barrier and the transmission coefficient 
depend on which RC is taken, but the final result that 
combines the two outcomes is independent of this choice. 

The RF method as proved its value for many systems, 
but also has its drawbacks. Although its result is inde- 
pendent of the chosen RC, its efficiency does and sen- 
sitively determines its success or failure. A non-suitable 
choice of RC can result in hysteresis effects in the free en- 
ergy calculation, which frustrates an accurate estimation 
of the barrier. Besides, even if an accurate value for the 
free energy barrier can be obtained, the corresponding 
transmission coefficient will be very small and its evalua- 
tion will require an extremely large number of pathways. 
In practice, it has been experienced that finding a good 
RC can be extremely difficult in high dimensional com- 
plex systems. Notable examples are chemical reactions 
in solution, where the reaction mechanism often depends 
on highly non-trivial solvent rearrangements. Also, com- 
puter simulations of nucleation processes use very compli- 
cated order parameters to distinguish between particles 
belonging to the liquid and solid phase. This makes it 
unfeasible to construct a single RC that accurately de- 
scribes the exact place of cross-over transitions. As re- 
sult, hysteresis effects and low transmission coefficients 
are almost unavoidable. 

The problem of finding suitable RCs, has urged the de- 
velopment of alternative methods. In 1998, Dellago et al. 
came up with such an alternative method that they called 
transition path sampling (TPS) 51 ' 52 i 53 ' 54 . This approach 
can be described as a MC sampling of MD pathways. 
Using a detailed balance technique, a set of trajectories 



can be collected that satisfy some predetermined criteria. 
For instance, one can constrain the start- and end-point 
of the path in such a way that each trajectory connects 
the reactant and product state. An important point is 
that this sampling of successful reactive events does not 
require a RC that captures the reactive mechanism, but 
only needs an order parameter that can distinguish be- 
tween reactant and product state. In addition to this, 
the first series of TPS paper o 51 ' 52 i 53 i 54 also provided a 
route to calculate reaction rates. However, this approach 
has seldom been used due to its high computational cost. 
Moreover, within the context of the reaction rate calcu- 
lation, it is not so obvious to state that the TPS order 
parameter is actually very different from a RC. In this 
approach, the end-point of the path is forced to progress 
in successive steps from reactant to product state. Hence, 
the TPS order parameter needs to describe the interme- 
diate states as well just as a RC in the standard methods. 

Luckily, the algorithmic procedure to calculate reac- 
tion rates using the same path sampling framework im- 
proved considerably when the transition interface sam- 
pling (TIS) technique was devised 3 -^. TIS uses a flexible 
pathlength which reduces the number of required MD 
steps significantly. Moreover, the TIS method also elim- 
inates the need of so-called MC shifting moves that re- 
quired a considerable percentage of the simulation time 
in the TPS scheme. In addition, one can show that the 
new mathematical formulation of the reaction rate is less 
sensitive to recrossing events which guarantees a faster 
convergence. 

While TPS imposes conditions to the start- and end- 
point of the path to be within certain intervals of the 
RC, TIS imposes an interface crossing condition. Except 
for the technique proposed in Ref. [4JJ , TIS needs a RC 
just like the original TPS scheme. The RC is required 
to define a set of interfaces between the stable reactant 
and product states. However, unlike the standard RF 
methods, the TIS efficiency is relatively insensitive to 
the choice of RC as was first proven in [55| . This point 
is a strong advantage in complex systems where a 'good 
RC can be extremely difficult to find. 



B. the TIS algorithm 

The TIS algorithm works as follows. First step is to 
define a RC and a set of related values Ao,Ai,...,A„ 
with Xi < At+i. The subsets of phase- or configuration 
points for which the RC is exactly equal to A; basically 
define multidimensional surfaces or interfaces which give 
the name to this method. These values/interfaces should 
obey the following requirements: if the RC is lower than 
Ao = Aa, the system should be in the reactant state A; if 
the RC is higher than A„ = \b the system should be in 
the product state B; n and the positions for the interfaces 
in between, A^ with 1 = 1, 2, . . . , n — 1, should be set to 
optimize the efficiency. Further, the surface A^ should 
be set in such a way that whenever a MD simulation is 
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Reject trial move. Keep the old path. Count it again 
and restart procedure using the same path. 





I. Pick randomly a timcslicc of the 
old path [Like the red dot in b)]. 



II. Add small Gaussian distributed 
random changes dp's to the momenta 
of all particles. [Like in c) or d)] 

III. Continue with a probability 
mm{l,ea;p[/3(KIN (o) - KIN*"')]} 
'as in d)]. Otherwise Reject [as in c)] 



IV. Take a random number a from 
an uniform distribution [0 : 1]. Define 
the maximum allowed pathlength 
jjnax _ intjXt°Yo;], Continue. 



V. Integrate the Newtonian equations 
of motion backward in time until 

^- reaching \ {) [sec g)], X i+1 [see e)] or 
until the pathlength exceeds L max [f)]. 
Continue in case g). Otherwise reject. 

VI. Integrate the equations of motion 
. foi"ward in time until reaching Ao 

[h) and j)]. Ai+i [k)j, or when the 
whole path exceeds L max [i)]. 

+ m ; 

VII. Reject in case i) or when the 
path has no crossing with A, [case h)]. 

^- Otherwise, accept the new path. 
Replace the old path by the new 
one [Say k)]. Update statistics. 





^ Accept the new path. 
Repeat procedure from step I using this new path. 



FIG. 1: Illustration of the shooting algorithm in TIS. The panels a)-l) depict trajectories/trajectory-segments on a free energy 
surface. The dashed horizontal lines are the TIS interfaces (n = 4 and i — 2 in this case). The algorithm requires an initial 
path a) to start the loop. The length of this particular path l/ - 1 is eighteen timeslices (endpoints are not included). At step 
I, a random point is picked from this old path and some small randomized changes are applied to the velocities of all the 
particles (II), followed by a Metropolis acceptance/rejection step (III). In c) the new velocities have resulted in a much larger 
kinetic energy KIN*-™' and therefore this trial move is most likely rejected. Step IV is requir ed to maintain detailed balance 
between pathways of different lengths. For example, if the random number generator assigns a — .59 then L max = 30 and we 
can reject when the path is unfinished, but already contains 31 timeslices as in panel f) and i). At V, the equations of motion 
are integrated backwards in time by a MD algorithm using the shooting point with reverse velocities as starting point. At VI) 
the equations of motion are integrated forward in time starting from the same shooting point (without reversed velocities). 
After a rejection the old path is kept and counted again. If accepted, the new path will automatically start at Xa and cross 
\i. The path can end at either Xa as in j) or at Xi+i as in k). The fraction of sampled pathways that end at Xi+i determines 
PA(A i+ i|Ai) 



released from within the reactant well, this surface should 
be frequently crossed. The TIS rate expression can then 
be formulated as 

k AB = f A V A (X B \X A ) (1) 



Here, f A is the escape flux through the first interface. 
In a long MD simulation, this simply corresponds to the 
number of detected crossings through the surface X A di- 
vided by the total simulation time (here, we assume that 
we will not observe a spontaneous transition to state B 
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during the simulation. For a more formal definition see 
[391]). The other term Va{^b\^a) is the overall crossing 
probability. This is the probability that whenever the 
system crosses A^, it will cross A_b before it crosses Xa 
again. As As is a surface at the other side of the bar- 
rier, this probability will be very small and can not be 
calculated directly. This probability can, however, be de- 
termined by a series of path sampling simulations using 
the following factorization: 

n-1 

Va{\b\Xa) = 7MA„|A ) = 1] v A(K+i\\i)- (2) 

2 = 

The terms Va{K+i \K) are history dependent conditional 
crossing probabilities which are much higher and can be 
computed. In words, P^(Ai+i|Ai) is the probability that 
Ai+i will be crossed before A,4 under the twofold condi- 
tion that the system is at the point to cross the interface 
A; in one timestep while Xa was more recently crossed 
than Xi in the past. It is due to this history dependence 
that Eq. $2$ is exact and should not be misinterpreted as 
a Markovian approximation. Va^X^i |Ai) is also equal to 
the number of all possible paths that start at Xa and end 
at Xi + i divided by the number of all possible paths that 
start at Xa, end at either A;+i or A ,4, and have at least 
one crossing with A^. Hence, this term can be calculated 
if we can generate the appropriate trajectories with their 
correct statistical weight. It is however not so obvious 
to generate these pathways especially when A^ is in the 
reaction barrier region. This difficulty can be overcome 
by a MC algorithm that employs a variation of the TPS 
shooting move. The algorithm is explained in Fig. [TJ 

The algorithm requires to have one path fulfilling the 
correct condition. That is starting at Xa and crossing 
Xi at least once before ending at either A^+i or A^. A 
crucial point is step II. After picking a random timeslice 
(a point that constitutes all the particle positions and 
momenta at a certain timestep along the path) , one adds 
random values to all the momenta. In practice, these ran- 
dom values are taken from a Gaussian distribution with 
a certain width a, that should be adapted to obtain the 
optimum efficiency. If a is small, the random momentum 
changes will be small as well and the new path will lie 
closely to the old one (if we assume deterministic dynam- 
ics). This small deviation results in a significant chance 
that the trial path will satisfy the required conditions as 
well which yields a good acceptance rate. However, a 
too small value of a will result in too strong correlations 
between the accepted moves (In the extreme case when 
(7 = 0, one regenerates exclusively the same path). Usu- 
ally, one tries several values for a in a series of short test 
simulations. It is generally assumed that the value that 
yields an acceptance rate of 50 % is close to an optimum 
value for a. If one wants to simulate at constant energy 
instead of constant temperature, step III can be replaced 
by a proper velocity rescaling procedure that, if needed, 
can also preserve linear and angular momentum 5 ^. 

Another important point is that, in order to enter the 



loop, one needs to have a single path that obeys the cor- 
rect requirements. This can already be quite difficult and 
several techniques to get such a first initial path have 
been suggested 5 -^,. However, in TIS these initial paths 
are generated automatically when the different types of 
simulations are consecutively performed (See Fig. [2J). 

First the MD simulation is performed to calculate Ja- 
Then, a series of path-sampling simulations follows to cal- 
culate V A (X 1 \Xo),VA(X2\X 1 ),...,V A (X n \X n - 1 ). When 
these simulations are performed in this order, each path- 
sampling simulation can obtain the necessary initial path 
from the previous simulation (See Fig. \2§- 

It is important to note that the final result, the reac- 
tion rate Uab, does not sensitively depend on the posi- 
tions of the outer interfaces A^ and As as long as they 
are reasonable. The number of interfaces n and their 
positions only influence the efficiency of the method. It 
was found that the total efficiency is optimized when for 
each path-simulation one out of five trajectories reaches 
the next interfac e) 41 ! 55 . Hence, using some initial trial 
simulations, one can adjust the number of interfaces and 
their position to satisfy this condition. The easiest way 
to achieve this is to use a slight variation of the algorithm 
that is shown in Fig. [TJ Instead of stopping the integra- 
tion when the trajectory crosses A^+i as in panel e) and 
k), one can continue the trajectory until it reaches Xa or 
A b ■ This algorithm only requires knowing the position of 
Xa, A_b, and A^. By examining the progress of the paths 
along the RC beyond Xi, one can define the next interface 
Ai + i exactly at the point where 80 % of the paths have 
returned to A^. 



C. analysis of the reaction mechanism 

When the complete series of simulations is finished, the 
reaction rate follows simply from Eqs. (|II2p . In addition 
to this, the ensemble of pathways can be analyzed which 
can yield valuable information about the reaction mech- 
anism. In this respect, the TIS path-ensembles might 
actually prove to be more useful than the ones obtained 
by the original TPS method. As each ensemble contains 
the correct ratio of paths progressing upto a certain level, 
but then either return or make a little step further, one 
can try to understand the characteristic differences be- 
tween the 'successful' and unsuccessful' pathways. In 
contrast, the TPS method aims to generate successful 
trajectories only. One of the properties that can improve 
understanding of mechanisms is the overall crossing prob- 
ability function. This represents a sort of path survival 
probability along the RC. This function equals 1 at Ao 
and Pa(Xb\Xa) at A^. In transition, this function is 
monotonically decreasing and terminates in a horizon- 
tal plateau when the barrier ridge is crossed completely. 
This function could be considered as a dynamical equiva- 
lent of the free energy profile along the RC. In Fig. [3] the 
overall crossing probability function is depicted together 
with the free energy profile obtained from a nucleation 
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process of Lennard- Jones particles. 



FIG. 2: This pictures illustrates some typical trajectories on 
a free energy surface for a series of TIS simulations. The 
glassy plates represent the TIS interfaces. The first type of 
simulation (top) , is a straightforward MD simulation which is 
required to calculate the flux /a through the first interface. 
The next step is a path-sampling simulation which generates 
pathways that start at A a and end at either A a or Ai. This 
simulation yields the result of 'Pa(AiIAo) and is illustrated 
in the middle panel. The initial path to start this simula- 
tion can be obtained by taking a trajectory segment from 
the initial MD simulation. The bottom panel shows the next 
path-sampling simulation to calculate 7- , a(A2|Ai). It gener- 
ates pathways that start at Xa and cross Ai at least once. 
Also here, the initial path can be obtained from the previ- 
ous simulation by taking one of the paths that successfully 
reached Ai. 



1 

0.1 

>, 

= 0.01 

\a 
ra 

-§ 0.001 

Q. 

P 1 e-04 

Tn 
to 

O 1e-05 




30 
25 
20 
15 

10^ 

5 



100 200 300 400 500 600 700 800 

largest solid cluster 

FIG. 3: The crossing probability function obtained by TIS 
and partial path TIS (PPTIS) 59 for nucleation process of LJ 
particles. The simulation data are obtained from [58]. The 
free energy profile, that was calculated simultaneously using 
the technique of [72j], is also shown. The simulated system 
contained 10648 particles in total. The RC was denned as the 
largest solid cluster in the system. More details can be found 
in [58l ]. The final crossing probability was Pa (A_b|Aa) = 
(7.8 ± 5.5) ■ I(T 7 (TIS) and (14.1 ± 8.7) • 1(T 7 (PPTIS). The 
crossing probability shows a plateau at 410.5 indicating that 
the reaction barrier ridge is crossed. The free energy barrier 
has a maximum at 243. After this point about 50 % (PPTIS) 
till 75 % (TIS) of the paths still fail to reach the reactant 
state. Although, the PPTIS and TIS results are within each 
others error-bars, it is likely that PPTIS overestimates the 
crossing probability due to the imprecise choice of RC 55 . 



The results were obtained from from Ref. [58|. In- 
terestingly, after the maximum of free energy barrier is 
crossed (cluster size 243), the majority of the trajectories 
(~ 75%) still fail to reach the reactant state. This effect 
might be partly due to diffusive motion, but is most likely 
an effect of an improperly chosen RC. This shows that 
the projection on a single RC using static free energy 
calculations can be misleading. Neither the height of the 
barrier nor the position of the TS dividing surface have 
to reflect the actual height and position of the reaction 
barrier. Indeed, Moroni et al. found that a 'good RC 
should at least incorporate one more important quantity 
which is the crystallinity of the cluster. Small clusters 
< 243 with a high crystallinity were found to grow fur- 
ther easily, while large clusters with less structure were 
unstable and broke-up into smaller pieces. 



D. recent and future developments and 
applications 

Besides an efficient algorithm for the calculation of re- 
action rates, the TIS method has provided a new math- 
ematical framework to describe rare events. Recent new 
simulation techniques have exploited this TIS theory. 
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For example, for diffusive barrier crossings, where transi- 
tion paths become very long, the partial path TIS (PP- 
TIS) method was devised 59 . Here, using the assump- 
tion of memory loss, much shorter paths are generated 
after which the overall crossing probability can be recon- 
structed by a recursive formulation. The forward flux 
sampling techniqu e) 60 ' 61 is basically the same as TIS, but 
the way to generate pathways is different. In this ap- 
proach, the endpoints of all the pathways successfully 
reaching the next interface are stored and starting from 
each point a set of new pathways is generated in the next 
simulation. The main advantage is that the FFS scheme 
supplies a route to handle non-equilibrium systems. A 
disadvantage is that this only works for stochastic dy- 
namics and will always yield much stronger correlations 
between the generated pathways and the different path 
ensembles even for the pure Brownian dynamics case. 
Another drawback of PPTIS and FFS methods is that 
they do not possess the same RC insensitivity as TIS 55 . 
If necessary, even a fully RC free approach is possible 
as was suggested in (4f[ using the TIS pathlength as a 
transition parameter. A nice feature of this approach 
is that it does not require to specify a specific product 
state. Combinations with Configurational Bias M C 41 ' 61 
and path swapping technique o 41 ' 62 may also yield promis- 
ing advances for the computational efficiency. 

The TIS and its variations have been applied to 
various systems ranging from simple test-system a 39 ' 59 , 
nucleatio n 58 ' 63 ' 64 , protein foldin g 65 ' 66 , biochemical net- 
works 69 , driven polymer translocation through pores 61 , 
micelle formation 67 , ab initio simulation of chemical re- 
actions^ - and DNA denaturation 6 ^. 

The TIS technique can open many possible avenues 
in the field of zeolite formation simulations at several 
stages of the process. For instance, the first elementary 
step to Si polymerization is the condensation reaction 
2Si(OH) 4 -> Si 2 0(OH) 6 +H 2 0. This has been studied by 
ab initio static analysis including implicit solvent 6 - 8 -. This 
study has revealed two possible reaction mechanisms. 
Such a system is small enough to be treated by ab initio 
MD 37 ' 38 including explicit solvent molecules. The TIS 
method could hence give valuable insight which reaction 
mechanism dominates when dynamics and explicit sol- 
vent is taken into account. Classical reactive forcefields 
and rare event methods, such as TIS and PPTIS, should 
make it possible to simulate the dynamics of Si polymer- 
ization at much lower temperatures than hitherto was 
possible 69 . This allows to study this process at condi- 
tions that are much closer to the experimental situation. 
Moreover, by a right construction of the interfaces, the 
TIS method allows to focus on reaction mechanisms and 
rates of some very specific polymerization reactions, for 
instance the formation of the Si-30/33 precursor—. 

With the development of lattice and KMC models, the 
study of the next stages of nucleation and zeolite growth 
also come into reach. The combination of KMC and path 
sampling is a promising yet unexplored territory. De- 
spite the enormous long simulation periods that can be 



achieved by KMC, the expectation time to form a criti- 
cal nucleus starting from a disorder solution is generally 
still out of reach. Therefore, most KMC studies have 
concentrated on growth rather than nucleation. Hence, 
the study of zeolite nucleation might benefit significantly 
using combined KMC and path sampling techniques. 



III. SUMMARY 

We have reviewed the TIS method, its variations and 
its possible applications for the theoretical study of zeo- 
lite synthesis. The TIS method is a an elegant approach 
circumventing the timescale problem not by speeding up 
the dynamics of the system itself, but by concentrating 
on the short time trajectories which are of interest with- 
out using any approximation. TIS allows to overcome 
reaction barriers by a sequence of simulation series. It 
is important to realize that the barrier crossing event is 
not enhanced due to some artifical force but only due 
to the MC acception/rejection steps that include the in- 
terface crossing condition. Hence, each trajectory in the 
TIS path ensembles satisfy the correct dynamics on the 
true potential energy surface. This makes the method 
fundamentally different from, for instance, the metady- 
namics 70 approach. The TIS method makes use of the 
fact that the time needed to actually cross the barrier, 
the transition time, is much shorter than relaxation time 
^ab ' w hich is the time wherein one can expect a reactive 
event from an arbitrary point within the reactant well. 

TIS can be combined with any type of dynamics such 
as ab initio MD, Langevin, pure Brownian motion, classi- 
cal MD and KMC. A requirement for application of this 
method is that the simulation of short trajectories can 
occur sufficiently fast. This limits the size of the systems 
which can be studied, ranging from several molecules for 
ab initio dynamics to several thousand molecules for MD, 
and even larger assemblies for KMC simulations. 

Still, substantial work has to be done before fully re- 
alistic modeling of zeolite synthesis is our reach. An 
important requirement is the development of more ac- 
curate reactive force fields that can describe chemical 
events within the environment of solvent and template 
molecules. Recently, a more systematic approach for this 
development was suggested^ 1 -. Even though, the lattice 
and KMC models are making substantial steps forward, 
inclusion of solvent effects in a lattice-type models has 
proven to be a difficult problem that has not yet been 
solved. To conclude, the simulation methods have made 
prodigious advancements in recent years and might ul- 
timately give answers to important questions regarding 
zeolite synthesis, that can not be unambiguously accessed 
by experimental techniques. The TIS methods can help 
in obtaining dynamical information for the crucial but 
rare reaction steps in the zeolite process. In the near fu- 
ture, we are going to explore the application of TIS for 
the study of zeolite genesis. 
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